## Code to create data for Figure 1a and Figure 1b
## Stability of Immigration Preferences
## Alexander Kustov, Dillon Laaker, and Cassidy Reller
remove(list = ls())
library("psych")
library("survey")
library("foreign")

setwd("/Users/dillonlaaker/Box Sync/Stability/data/datasets")
ines <- read.dta("ines_data.dta")
bes <- read.dta("bes_data.dta")
taps <- read.dta("taps_data.dta")
ncp <- read.dta("ncp_data.dta")
shp <- read.dta("shp_data.dta")
liss <- read.dta("liss_data.dta")
ess <- read.dta("ess.dta")
anes <- read.dta("anes.dta")
gles <- read.dta("gles_data.dta")
cces <- read.dta("cces_data.dta")
vsg <- read.dta("vsg_data.dta")

##ESS
attach(ess)
ess <- data.frame(cntry, imsmetn, imdfetn, impcntr, imbgeco, imueclt, imwbcnt, essround, dweight, pweight)
detach(ess)
ess$weight <- ess$dweight*ess$pweight

ess$imm1 <- (ess$imsmetn-1)/3
ess$imm2 <- (ess$imdfetn-1)/3
ess$imm3 <- (ess$impcntr-1)/3
ess$imm4 <- (ess$imbgeco)/10*(-1)+1
ess$imm5 <- (ess$imueclt)/10*(-1)+1
ess$imm6 <- (ess$imwbcnt)/10*(-1)+1
ess$index <- (ess$imm1+ess$imm2+ess$imm3+ess$imm4+ess$imm5+ess$imm6)/6
ess$country[ess$cntry=="AT"] <- "Austria"
ess$country[ess$cntry=="BE"] <- "Belgium"
ess$country[ess$cntry=="CH"] <- "Switzerland"
ess$country[ess$cntry=="DE"] <- "Germany"
ess$country[ess$cntry=="DK"] <- "Denmark"
ess$country[ess$cntry=="ES"] <- "Spain"
ess$country[ess$cntry=="FI"] <- "Finland"
ess$country[ess$cntry=="FR"] <- "France"
ess$country[ess$cntry=="GB"] <- "Britain"
ess$country[ess$cntry=="GR"] <- "Greece"
ess$country[ess$cntry=="IE"] <- "Ireland"
ess$country[ess$cntry=="IL"] <- "Israel"
ess$country[ess$cntry=="IS"] <- "Iceland"
ess$country[ess$cntry=="IT"] <- "Italy"
ess$country[ess$cntry=="NL"] <- "Netherlands"
ess$country[ess$cntry=="NO"] <- "Norway"
ess$country[ess$cntry=="PL"] <- "Poland"
ess$country[ess$cntry=="PT"] <- "Portugal"
ess$country[ess$cntry=="RU"] <- "Russia"
ess$country[ess$cntry=="SE"] <- "Sweden"
ess$date[ess$essround==1] <- "01/01/2002"
ess$date[ess$essround==2] <- "01/01/2004"
ess$date[ess$essround==3] <- "01/01/2006"
ess$date[ess$essround==4] <- "01/01/2008"
ess$date[ess$essround==5] <- "01/01/2010"
ess$date[ess$essround==6] <- "01/01/2012"
ess$date[ess$essround==7] <- "01/01/2014"
ess$date <- as.Date(ess$date, format="%m/%d/%Y")
dsgness <- svydesign(ids=~1, data=ess, weights=~dweight)
y <- svyby(~index, ~date+country, dsgness, svymean, na.rm=TRUE)
z <- data.frame(y)
z$survey <- "ESS"
z$imm <- z$index
attach(z)
essmeans <- data.frame(survey, country, date, imm)
detach(z)
attach(essmeans)
essmeans1 <- essmeans[country=="Netherlands" | country=="Austria" |  country=="Germany" | country=="Russia", ]
essmeans2 <- essmeans[country=="Belgium" | country=="Italy" | country=="Finland" | country=="Iceland", ]
essmeans3 <- essmeans[country=="Greece" | country=="Britain" |  country=="Switzerland" | country=="Spain", ]
essmeans4 <- essmeans[country=="Norway" | country=="France" | country=="Sweden" | country=="Portugal", ]
essmeans5 <- essmeans[country=="Ireland" | country=="Denmark" |  country=="Poland" | country=="Israel", ]
detach(essmeans)
dsgness <- svydesign(ids=~1, data=ess, weights=~weight)
y  <- svyby(~index, ~date, dsgness, svymean, na.rm=TRUE)
z  <-  data.frame(y)
z$survey <- "ESS"
z$country <- "ESS"
z$imm <- z$index
attach(z)
essmeans <- data.frame(survey, country, date, imm)
detach(z)

data <- data.frame(ess$imm1, ess$imm2, ess$imm3, ess$imm4, ess$imm5, ess$imm6)
alpha(data) ## Cronbach's Alpha for ESS questions - measures internal consistency of data

##ANES
attach(anes)
anes <- data.frame(vcf0009z, vcf0879, vcf0004)
detach(anes)
anes$weight <- anes$vcf0009z
anes$imm <- anes$vcf0879
anes$imm[anes$vcf0879==8] <- NA
anes$imm[anes$vcf0879==9] <- NA
anes$imm <- (anes$imm-1)/4
anes$year[anes$vcf0004==1992] <- 1992
anes$year[anes$vcf0004==1994] <- 1994
anes$year[anes$vcf0004==1996] <- 1996
anes$year[anes$vcf0004==1998] <- 1998
anes$year[anes$vcf0004==2004] <- 2004
anes$year[anes$vcf0004==2008] <- 2008
anes$year[anes$vcf0004==2012] <- 2012
dsgnanes <- svydesign(ids=~1, data=anes, weights=~weight)
y  <- svyby(~imm, ~year, dsgnanes, svymean, na.rm=TRUE)
z  <-  data.frame(y)
z$country <- "US"
z$survey  <-  "ANES"
z$date <- c("01/01/1992", "01/01/1994", "01/01/1996", "01/01/1998", "01/01/2004", "01/01/2008", "01/01/2012")
z$date <- as.Date(z$date, format="%m/%d/%Y")
attach(z)
anesmeans  <-  data.frame(survey, country, date, imm)
detach(z)

##INES
inesmean2002 <- (mean(ines$imm2, na.rm=TRUE)-1)/6
inesmean2003 <- (mean(ines$imm3, na.rm=TRUE)-1)/6
inesmean2004 <- (mean(ines$imm4, na.rm=TRUE)-1)/6
inesmean2006 <- (mean(ines$imm6, na.rm=TRUE)-1)/6
inesmean2007 <- (mean(ines$imm7, na.rm=TRUE)-1)/6
imm <- c(inesmean2002, inesmean2003, inesmean2004, inesmean2006, inesmean2007)
date <- c("01/01/2002", "01/01/2003", "01/01/2004", "01/01/2006", "01/01/2007")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(5, "INES")
country <- replicate(5, "ines")
inesmeans <- data.frame(survey, country, date, imm)

inesmeanall2002 <- (mean(ines$imm2[ines$all==1], na.rm=TRUE)-1)/6
inesmeanall2003 <- (mean(ines$imm3[ines$all==1], na.rm=TRUE)-1)/6
inesmeanall2004 <- (mean(ines$imm4[ines$all==1], na.rm=TRUE)-1)/6
inesmeanall2006 <- (mean(ines$imm6[ines$all==1], na.rm=TRUE)-1)/6
inesmeanall2007 <- (mean(ines$imm7[ines$all==1], na.rm=TRUE)-1)/6
imm <- c(inesmeanall2002, inesmeanall2003, inesmeanall2004, inesmeanall2006, inesmeanall2007)
date <- c("01/01/2002", "01/01/2003", "01/01/2004", "01/01/2006", "01/01/2007")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(5, "INES - All Waves")
country <- replicate(5, "ines")
inesmeansall <- data.frame(survey, country, date, imm)

##LISS
lissmeans2008 <- (mean(liss$immigration_8, na.rm=TRUE)-1)/4
lissmeans2009 <- (mean(liss$immigration_9, na.rm=TRUE)-1)/4
lissmeans2010 <- (mean(liss$immigration_10, na.rm=TRUE)-1)/4
lissmeans2011 <- (mean(liss$immigration_11, na.rm=TRUE)-1)/4
lissmeans2012 <- (mean(liss$immigration_12, na.rm=TRUE)-1)/4
lissmeans2013 <- (mean(liss$immigration_13, na.rm=TRUE)-1)/4
lissmeans2014 <- (mean(liss$immigration_14, na.rm=TRUE)-1)/4
lissmeans2016 <- (mean(liss$immigration_16, na.rm=TRUE)-1)/4
lissmeans2017 <- (mean(liss$immigration_17, na.rm=TRUE)-1)/4
imm <- c(lissmeans2008, lissmeans2009, lissmeans2010, lissmeans2011, lissmeans2012, lissmeans2013, lissmeans2014, lissmeans2016, lissmeans2017)
date <- c("01/01/2008", "01/01/2009", "01/01/2010", "01/01/2011", "01/01/2012", "01/01/2013", "01/01/2014", "01/01/2016", "01/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(9, "LISS")
country <- replicate(9, "Netherlands")
lissmeans <- data.frame(survey, country, date, imm)

lissmeansall2008 <- (mean(liss$immigration_8[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2009 <- (mean(liss$immigration_9[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2010 <- (mean(liss$immigration_10[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2011 <- (mean(liss$immigration_11[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2012 <- (mean(liss$immigration_12[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2013 <- (mean(liss$immigration_13[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2014 <- (mean(liss$immigration_14[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2016 <- (mean(liss$immigration_16[liss$all==1], na.rm=TRUE)-1)/4
lissmeansall2017 <- (mean(liss$immigration_17[liss$all==1], na.rm=TRUE)-1)/4
imm <- c(lissmeansall2008, lissmeansall2009, lissmeansall2010, lissmeansall2011, lissmeansall2012, lissmeansall2013, lissmeansall2014, lissmeansall2016, lissmeansall2017)
date <- c("01/01/2008", "01/01/2009", "01/01/2010", "01/01/2011", "01/01/2012", "01/01/2013", "01/01/2014", "01/01/2016", "01/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(9, "LISS - All Waves")
country <- replicate(9, "Netherlands")
lissmeansall <- data.frame(survey, country, date, imm)


##SHP
shpmeans1999 <- (mean(shp$imm99, na.rm=TRUE)-1)/2
shpmeans2000 <- (mean(shp$imm0, na.rm=TRUE)-1)/2
shpmeans2001 <- (mean(shp$imm1, na.rm=TRUE)-1)/2
shpmeans2002 <- (mean(shp$imm2, na.rm=TRUE)-1)/2
shpmeans2003 <- (mean(shp$imm3, na.rm=TRUE)-1)/2
shpmeans2004 <- (mean(shp$imm4, na.rm=TRUE)-1)/2
shpmeans2005 <- (mean(shp$imm5, na.rm=TRUE)-1)/2
shpmeans2006 <- (mean(shp$imm6, na.rm=TRUE)-1)/2
shpmeans2007 <- (mean(shp$imm7, na.rm=TRUE)-1)/2
shpmeans2008 <- (mean(shp$imm8, na.rm=TRUE)-1)/2
shpmeans2009 <- (mean(shp$imm9, na.rm=TRUE)-1)/2
shpmeans2011 <- (mean(shp$imm11, na.rm=TRUE)-1)/2
imm <- c(shpmeans1999, shpmeans2000, shpmeans2001, shpmeans2002, shpmeans2003, shpmeans2004, shpmeans2005, shpmeans2006, shpmeans2007, shpmeans2008, shpmeans2009, shpmeans2011)
date <- c("01/01/1999", "01/01/2000", "01/01/2001", "01/01/2002", "01/01/2003", "01/01/2004", "01/01/2005", "01/01/2006", "01/01/2007", "01/01/2008", "01/01/2009", "01/01/2011")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(12, "SHP")
country <- replicate(12, "Switzerland")
shpmeans <- data.frame(survey, country, date, imm)

shpmeansall1999 <- (mean(shp$imm99[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2000 <- (mean(shp$imm0[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2001 <- (mean(shp$imm1[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2002 <- (mean(shp$imm2[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2003 <- (mean(shp$imm3[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2004 <- (mean(shp$imm4[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2005 <- (mean(shp$imm5[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2006 <- (mean(shp$imm6[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2007 <- (mean(shp$imm7[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2008 <- (mean(shp$imm8[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2009 <- (mean(shp$imm9[shp$all==1], na.rm=TRUE)-1)/2
shpmeansall2011 <- (mean(shp$imm11[shp$all==1], na.rm=TRUE)-1)/2
imm <- c(shpmeansall1999, shpmeansall2000, shpmeansall2001, shpmeansall2002, shpmeansall2003, shpmeansall2004, shpmeansall2005, shpmeansall2006, shpmeansall2007, shpmeansall2008, shpmeansall2009, shpmeansall2011)
date <- c("01/01/1999", "01/01/2000", "01/01/2001", "01/01/2002", "01/01/2003", "01/01/2004", "01/01/2005", "01/01/2006", "01/01/2007", "01/01/2008", "01/01/2009", "01/01/2011")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(12, "SHP - All Waves")
country <- replicate(12, "Switzerland")
shpmeansall <- data.frame(survey, country, date, imm)

besmeans1 <- (mean(bes$imm1, na.rm=TRUE)-1)/6
besmeans2 <- (mean(bes$imm2, na.rm=TRUE)-1)/6
besmeans3 <- (mean(bes$imm3, na.rm=TRUE)-1)/6
besmeans4 <- (mean(bes$imm4, na.rm=TRUE)-1)/6
besmeans7 <- (mean(bes$imm7, na.rm=TRUE)-1)/6
besmeans8 <- (mean(bes$imm8, na.rm=TRUE)-1)/6
besmeans10 <- (mean(bes$imm10, na.rm=TRUE)-1)/6
besmeans11 <- (mean(bes$imm11, na.rm=TRUE)-1)/6
imm <- c(besmeans1, besmeans2, besmeans3, besmeans4, besmeans7, besmeans8, besmeans10, besmeans11)
date <- c("02/01/2014", "05/01/2014", "09/01/2014", "03/01/2015", "04/01/2016", "05/01/2016", "11/01/2016", "04/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(8, "BES")
country <- replicate(8, "Britain")
besmeans <- data.frame(survey, country, date, imm)

besmeansall1 <- (mean(bes$imm1[bes$all==1], na.rm=TRUE)-1)/6
besmeansall2 <- (mean(bes$imm2[bes$all==1], na.rm=TRUE)-1)/6
besmeansall3 <- (mean(bes$imm3[bes$all==1], na.rm=TRUE)-1)/6
besmeansall4 <- (mean(bes$imm4[bes$all==1], na.rm=TRUE)-1)/6
besmeansall7 <- (mean(bes$imm7[bes$all==1], na.rm=TRUE)-1)/6
besmeansall8 <- (mean(bes$imm8[bes$all==1], na.rm=TRUE)-1)/6
besmeansall10 <- (mean(bes$imm10[bes$all==1], na.rm=TRUE)-1)/6
besmeansall11 <- (mean(bes$imm11[bes$all==1], na.rm=TRUE)-1)/6
imm <- c(besmeansall1, besmeansall2, besmeansall3, besmeansall4, besmeansall7, besmeansall8, besmeansall10, besmeansall11)
date <- c("02/01/2014", "05/01/2014", "09/01/2014", "03/01/2015", "04/01/2016", "05/01/2016", "11/01/2016", "04/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(8, "BES - All Waves")
country <- replicate(8, "Britain")
besmeansall <- data.frame(survey, country, date, imm)

ncpmeans3 <- (mean(ncp$immigration3, na.rm=TRUE)-1)/6
ncpmeans4 <- (mean(ncp$immigration4, na.rm=TRUE)-1)/6
ncpmeans5 <- (mean(ncp$immigration5, na.rm=TRUE)-1)/6
ncpmeans6 <- (mean(ncp$immigration6, na.rm=TRUE)-1)/6
ncpmeans7 <- (mean(ncp$immigration7, na.rm=TRUE)-1)/6
ncpmeans8 <- (mean(ncp$immigration8, na.rm=TRUE)-1)/6
imm <- c(ncpmeans3, ncpmeans4, ncpmeans5, ncpmeans6, ncpmeans7, ncpmeans8)
date <- c("10/01/2014", "03/01/2015", "11/01/2015", "03/01/2016", "11/01/2016", "03/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(6, "NCP")
country <- replicate(6, "ncp")
ncpmeans <- data.frame(survey, country, date, imm)

ncpmeansall3 <- (mean(ncp$immigration3[ncp$all==1], na.rm=TRUE)-1)/6
ncpmeansall4 <- (mean(ncp$immigration4[ncp$all==1], na.rm=TRUE)-1)/6
ncpmeansall5 <- (mean(ncp$immigration5[ncp$all==1], na.rm=TRUE)-1)/6
ncpmeansall6 <- (mean(ncp$immigration6[ncp$all==1], na.rm=TRUE)-1)/6
ncpmeansall7 <- (mean(ncp$immigration7[ncp$all==1], na.rm=TRUE)-1)/6
ncpmeansall8 <- (mean(ncp$immigration8[ncp$all==1], na.rm=TRUE)-1)/6
imm <- c(ncpmeansall3, ncpmeansall4, ncpmeansall5, ncpmeansall6, ncpmeansall7, ncpmeansall8)
date <- c("10/01/2014", "03/01/2015", "11/01/2015", "03/01/2016", "11/01/2016", "03/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(6, "NCP - All Waves")
country <- replicate(6, "ncp")
ncpmeansall <- data.frame(survey, country, date, imm)

tapsmeans8 <- mean(taps$imm8, na.rm=TRUE)
tapsmeans15 <- mean(taps$imm15, na.rm=TRUE)
tapsmeans22 <- mean(taps$imm22, na.rm=TRUE)
tapsmeans27 <- mean(taps$imm27, na.rm=TRUE)
tapsmeans34 <- mean(taps$imm34, na.rm=TRUE)
tapsmeans37 <- mean(taps$imm37, na.rm=TRUE)
tapsmeans46 <- mean(taps$imm46, na.rm=TRUE)
tapsmeans49 <- mean(taps$imm49, na.rm=TRUE)
tapsmeans51 <- mean(taps$imm51, na.rm=TRUE)
tapsmeans54 <- mean(taps$imm54, na.rm=TRUE)
tapsmeans56 <- mean(taps$imm56, na.rm=TRUE)
imm <- c(tapsmeans8, tapsmeans15, tapsmeans22, tapsmeans27, tapsmeans34, tapsmeans37, tapsmeans46, tapsmeans49, tapsmeans51, tapsmeans54, tapsmeans56)
date <- c("07/01/2012", "02/01/2013", "09/01/2013", "02/01/2014", "09/01/2014", "12/01/2014", "09/01/2015", "12/01/2015", "02/01/2016", "5/01/2016", "07/01/2016")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(11, "TAPS")
country <- replicate(11, "U.S.")
tapsmeans <- data.frame(survey, country, date, imm)

tapsmeansall8 <- mean(taps$imm8[taps$all==1], na.rm=TRUE)
tapsmeansall15 <- mean(taps$imm15[taps$all==1], na.rm=TRUE)
tapsmeansall22 <- mean(taps$imm22[taps$all==1], na.rm=TRUE)
tapsmeansall27 <- mean(taps$imm27[taps$all==1], na.rm=TRUE)
tapsmeansall34 <- mean(taps$imm34[taps$all==1], na.rm=TRUE)
tapsmeansall37 <- mean(taps$imm37[taps$all==1], na.rm=TRUE)
tapsmeansall46 <- mean(taps$imm46[taps$all==1], na.rm=TRUE)
tapsmeansall49 <- mean(taps$imm49[taps$all==1], na.rm=TRUE)
tapsmeansall51 <- mean(taps$imm51[taps$all==1], na.rm=TRUE)
tapsmeansall54 <- mean(taps$imm54[taps$all==1], na.rm=TRUE)
tapsmeansall56 <- mean(taps$imm56[taps$all==1], na.rm=TRUE)
imm <- c(tapsmeansall8, tapsmeansall15, tapsmeansall22, tapsmeansall27, tapsmeansall34, tapsmeansall37, tapsmeansall46, tapsmeansall49, tapsmeansall51, tapsmeansall54, tapsmeansall56)
date <- c("07/01/2012", "02/01/2013", "09/01/2013", "02/01/2014", "09/01/2014", "12/01/2014", "09/01/2015", "12/01/2015", "02/01/2016", "5/01/2016", "07/01/2016")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(11, "TAPS - All Waves")
country <- replicate(11, "U.S.")
tapsmeansall <- data.frame(survey, country, date, imm)

glesmeans1 <- (mean(gles$imm1, na.rm=TRUE) - 1)/6
glesmeans2 <- (mean(gles$imm2, na.rm=TRUE) - 1)/6
glesmeans4 <- (mean(gles$imm4, na.rm=TRUE) - 1)/6
glesmeans6 <- (mean(gles$imm6, na.rm=TRUE) - 1)/6
glesmeans10 <- (mean(gles$imm10, na.rm=TRUE) - 1)/6
glesmeans11 <- (mean(gles$imm11, na.rm=TRUE) - 1)/6
glesmeans12 <- (mean(gles$imm12, na.rm=TRUE) - 1)/6
glesmeans13 <- (mean(gles$imm13, na.rm=TRUE) - 1)/6
glesmeans15 <- (mean(gles$imm15, na.rm=TRUE) - 1)/6
glesmeans17 <- (mean(gles$imm17, na.rm=TRUE) - 1)/6
imm <- c(glesmeans1, glesmeans2, glesmeans4, glesmeans6, glesmeans10, glesmeans11, glesmeans12, glesmeans13, glesmeans15, glesmeans17)
date <- c("06/01/2013", "07/01/2013", "08/01/2013", "09/01/2013", "06/01/2016", "02/01/2017", "05/01/2017", "07/01/2017", "09/01/2017", "10/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(10, "GLES")
country <- replicate(10, "Germany")
glesmeans <- data.frame(survey, country, date, imm)

glesmeansall1 <- (mean(gles$imm1[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall2 <- (mean(gles$imm2[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall4 <- (mean(gles$imm4[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall6 <- (mean(gles$imm6[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall10 <- (mean(gles$imm10[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall11 <- (mean(gles$imm11[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall12 <- (mean(gles$imm12[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall13 <- (mean(gles$imm13[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall15 <- (mean(gles$imm15[gles$all==1], na.rm=TRUE) - 1)/6
glesmeansall17 <- (mean(gles$imm17[gles$all==1], na.rm=TRUE) - 1)/6
imm <- c(glesmeansall1, glesmeansall2, glesmeansall4, glesmeansall6, glesmeansall10, glesmeansall11, glesmeansall12, glesmeansall13, glesmeansall15, glesmeansall17)
date <- c("06/01/2013", "07/01/2013", "08/01/2013", "09/01/2013", "06/01/2016", "02/01/2017", "05/01/2017", "07/01/2017", "09/01/2017", "10/01/2017")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(10, "GLES - All Waves")
country <- replicate(10, "Germany")
glesmeansall <- data.frame(survey, country, date, imm)

ccesmeans10 <- (mean(cces$immigration_10, na.rm=TRUE))
ccesmeans12 <- (mean(cces$immigration_12, na.rm=TRUE))
ccesmeans14 <- (mean(cces$immigration_14, na.rm=TRUE))
imm <- c(ccesmeans10, ccesmeans12, ccesmeans14)
date <- c("01/01/2010", "01/01/2012", "01/01/2014")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(3, "CCES")
country <- replicate(3, "US")
ccesmeans <- data.frame(survey, country, date, imm)

ccesmeansall10 <- (mean(cces$immigration_10[cces$all==1], na.rm=TRUE))
ccesmeansall12 <- (mean(cces$immigration_12[cces$all==1], na.rm=TRUE))
ccesmeansall14 <- (mean(cces$immigration_14[cces$all==1], na.rm=TRUE))
imm <- c(ccesmeansall10, ccesmeansall12, ccesmeansall14)
date <- c("01/01/2010", "01/01/2012", "01/01/2014")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(3, "CCES - All Waves")
country <- replicate(3, "US")
ccesmeansall <- data.frame(survey, country, date, imm)

vsgmeans11 <- (mean(vsg$immigration_11, na.rm=TRUE))
vsgmeans16 <- (mean(vsg$immigration_16, na.rm=TRUE))
vsgmeans18 <- (mean(vsg$immigration_18, na.rm=TRUE))
imm <- c(vsgmeans11, vsgmeans16, vsgmeans18)
date <- c("01/01/2011", "01/01/2016", "01/01/2018")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(3, "VSG")
country <- replicate(3, "US")
vsgmeans <- data.frame(survey, country, date, imm)

vsgmeansall11 <- (mean(vsg$immigration_11[vsg$all==1], na.rm=TRUE))
vsgmeansall16 <- (mean(vsg$immigration_16[vsg$all==1], na.rm=TRUE))
vsgmeansall18 <- (mean(vsg$immigration_18[vsg$all==1], na.rm=TRUE))
imm <- c(vsgmeansall11, vsgmeansall16, vsgmeansall18)
date <- c("01/01/2011", "01/01/2016", "01/01/2018")
date <- as.Date(date, format="%m/%d/%Y")
survey <- replicate(3, "VSG - All Waves")
country <- replicate(3, "US")
vsgmeansall <- data.frame(survey, country, date, imm)

figureA1a <- rbind(anesmeans, inesmeans, inesmeansall, lissmeans, lissmeansall, shpmeans, shpmeansall, vsgmeans, vsgmeansall, essmeans)
figureA1b <- rbind(besmeans, besmeansall, ccesmeans, ccesmeansall, ncpmeans, ncpmeansall, tapsmeans, tapsmeansall, glesmeans, glesmeansall)
figureA3a <- essmeans1
figureA3b <- essmeans2
figureA3c <- essmeans3
figureA3d <- essmeans4
figureA3e <- essmeans5

setwd("/Users/dillonlaaker/Box Sync/Stability/data/figures/figureA1")
write.csv(figureA1a, file="figureA1a.csv")
write.csv(figureA1b, file="figureA1b.csv")

remove(list = ls())

setwd("/Users/dillonlaaker/Box Sync/Stability/data/figures/figureA1")
figureA1a <- read.csv("figureA1a.csv")
figureA1b <- read.csv("figureA1b.csv")
figureA1a$date <- as.Date(figureA1a$date, format = "%Y-%m-%d")
figureA1b$date <- as.Date(figureA1b$date, format = "%Y-%m-%d")

setwd("/Users/dillonlaaker/Box Sync/Stability/Draft")

dates <- c("01/01/1990", "01/01/1991", "01/01/1992", "01/01/1993", "01/01/1994", "01/01/1995", "01/01/1996", "01/01/1997", "01/01/1998", "01/01/1999", "01/01/2000", "01/01/2001", "01/01/2002", "01/01/2003", "01/01/2004", "01/01/2005", "01/01/2006", "01/01/2007", "01/01/2008", "01/01/2009", "01/01/2010", "01/01/2011", "01/01/2012", "01/01/2013", "01/01/2014", "01/01/2015", "01/01/2016", "01/01/2017", "01/01/2018")
dates <- as.Date(dates, format = "%m/%d/%Y")
xlimd <- c("01/01/1990", "03/01/2018") 
xlimd <- as.Date(xlimd, format = "%m/%d/%Y")
pdf("figureA1a.pdf", height=6, width=12)
par(mgp=c(3,1,0))
par(mar=c(5,5,1,1.5))
plot(figureA1a$date[figureA1a$survey=="ANES"], figureA1a$imm[figureA1a$survey=="ANES"], type="o", lty=1, pch=10, bg="white", ylim=c(0,1), yaxt="n", xaxt="n", bty="l", xaxs="i", yaxs="i", xlim=xlimd, ylab="", xlab="")
mtext(side=2, text="Immigration Score (0-1)", line=3, cex=1.25)
mtext(side=1, text="Year", line=3, cex=1.25)
axis(side=2, at=c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1), labels=c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1), lty=1, las=2)
axis.Date(side=1, at=dates, label=c("1990", "", "", "", "", "1995", "", "", "", "", "2000", "", "", "", "", "2005", "", "", "", "", "2010", "", "", "", "", "2015", "", "", ""), lty=1, las=0)
lines(figureA1a$date[figureA1a$survey=="LISS"], figureA1a$imm[figureA1a$survey=="LISS"], type="o", lty=1, pch=22, bg="white")
lines(figureA1a$date[figureA1a$survey=="LISS - All Waves"], figureA1a$imm[figureA1a$survey=="LISS - All Waves"], type="o", lty=2, pch=22, bg="white")
lines(figureA1a$date[figureA1a$survey=="VSG"], figureA1a$imm[figureA1a$survey=="VSG"], type="o", lty=1, pch=8, bg="white")
lines(figureA1a$date[figureA1a$survey=="VSG - All Waves"], figureA1a$imm[figureA1a$survey=="VSG - All Waves"], type="o", lty=2, pch=8, bg="white")
lines(figureA1a$date[figureA1a$survey=="SHP"], figureA1a$imm[figureA1a$survey=="SHP"], type="o", lty=1, pch=9, bg="white")
lines(figureA1a$date[figureA1a$survey=="SHP - All Waves"], figureA1a$imm[figureA1a$survey=="SHP - All Waves"], type="o", lty=2, pch=9, bg="white")
lines(figureA1a$date[figureA1a$survey=="INES"], figureA1a$imm[figureA1a$survey=="INES"], type="o", lty=1, pch=24, bg="white")
lines(figureA1a$date[figureA1a$survey=="INES - All Waves"], figureA1a$imm[figureA1a$survey=="INES - All Waves"], type="o", lty=2, pch=24, bg="white")
lines(figureA1a$date[figureA1a$survey=="ESS"], figureA1a$imm[figureA1a$survey=="ESS"], type="o", lty=1, pch=7, bg="white")
legend("bottomright", legend=c("ANES", "ESS", "INES", "LISS", "SHP", "VSG"), pch=c(10, 7, 24, 22, 9, 8), bty="n", pt.cex=1.25, cex=1.25)
dev.off()

dates <- c("12/01/2009", "06/01/2010", "12/01/2010", "06/01/2011", "12/01/2011", "06/01/2012", "12/01/2012", "06/01/2013", "12/01/2013", "06/01/2014", "12/01/2014", "06/01/2015", "12/01/2015", "06/01/2016", "12/01/2016", "06/01/2017", "12/01/2017")
dates <- as.Date(dates, format = "%m/%d/%Y")
xlimd <- c("12/01/2009", "12/01/2017") 
xlimd <- as.Date(xlimd, format = "%m/%d/%Y")
pdf("figureA1b.pdf", height=6, width=12)
par(mgp=c(3,1,0))
par(mar=c(5,5,1,1.5))
plot(figureA1b$date[figureA1b$survey=="BES"], figureA1b$imm[figureA1b$survey=="BES"], type="o", lty=1, pch=23, bg="white", ylim=c(0,1), yaxt="n", xaxt="n", bty="l", xaxs="i", yaxs="i", xlim=xlimd, ylab="", xlab="")
mtext(side=2, text="Immigration Score (0-1)", line=3, cex=1.25)
mtext(side=1, text="Month/Year", line=3, cex=1.25)
axis(side=2, at=c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1), labels=c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1), lty=1, las=2)
axis.Date(side=1, at=dates, label=c("12/2009", "", "12/2010", "", "12/2011", "", "12/2012", "", "12/2013", "", "12/2014", "", "12/2015", "", "12/2016", "", "12/2017"), lty=1, las=0)
lines(figureA1b$date[figureA1b$survey=="BES - All Waves"], figureA1b$imm[figureA1b$survey=="BES - All Waves"], type="o", lty=2, pch=23, bg="white")
lines(figureA1b$date[figureA1b$survey=="TAPS"], figureA1b$imm[figureA1b$survey=="TAPS"], type="o", lty=1, pch=4, bg="white")
lines(figureA1b$date[figureA1b$survey=="TAPS - All Waves"], figureA1b$imm[figureA1b$survey=="TAPS - All Waves"], type="o", lty=2, pch=4, bg="white")
lines(figureA1b$date[figureA1b$survey=="NCP"], figureA1b$imm[figureA1b$survey=="NCP"], type="o", lty=1, pch=21, bg="white")
lines(figureA1b$date[figureA1b$survey=="NCP - All Waves"], figureA1b$imm[figureA1b$survey=="NCP - All Waves"], type="o", lty=2, pch=21, bg="white")
lines(figureA1b$date[figureA1b$survey=="GLES"], figureA1b$imm[figureA1b$survey=="GLES"], type="o", lty=1, pch=25, bg="white")
lines(figureA1b$date[figureA1b$survey=="GLES - All Waves"], figureA1b$imm[figureA1b$survey=="GLES - All Waves"], type="o", lty=2, pch=25, bg="white")
lines(figureA1b$date[figureA1b$survey=="CCES"], figureA1b$imm[figureA1b$survey=="CCES"], type="o", lty=1, pch=11, bg="white")
lines(figureA1b$date[figureA1b$survey=="CCES - All Waves"], figureA1b$imm[figureA1b$survey=="CCES - All Waves"], type="o", lty=2, pch=11, bg="white")
legend("bottomright", legend=c("BES", "CCES", "GLES", "NCP", "TAPS"), pch=c(23, 11, 25, 21, 4), bty="n", pt.cex=1.25, cex=1.25)
dev.off()